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Abstract 

In the present paper we investigate the influence of inhomogeneities in the 
dynamics and stability of DNA open states, modeled as propagating solitons 
in the spirit of a Generalized Yakushevish Model. It is a direct consequence 
of our model that there exists a critical distance between the soliton's center 
of mass and the inhomogeneity at which the interaction between them can 
change the stability of the open state. Furtherly from this results was derived 
a renormalized potential function. 
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I. INTRODUCTION 



The origin and dynamics of open states in DNA have been the subject of intense study 
in the last three decades [|l|,0,[|f|,|5| . To date some nonlinear models have been proposed 
0JlJ3)@)Bll3S0 that addresses this issue, which is somehow related to the transcrip- 
tion or replication processes in this molecule, in terms of solitons and their interactions. 
The torsional model proposed by Yakushevich |5|J^,[^,^iP, |l0iPi^2ip!3| is one of such mod- 
els, and in this paper we generalize it to take into account the effects of inhomogeneity, 
which take place in real DNA molecules [|TJ],|rj|.This problem has been considered before 
|5|j6||^,3,Pi ro|lll|l2|,|13 , |14 , l5l, 16 1 but in the framework of a simpler version of the model; we 
are not aware that the effect of the inhomogeneities have been studied previously within the 
full Yakushevich model. 

This paper is organized as follows. In section 2 a description of our generalized model and 
the fundamental equations are presented. In section 3 an analysis of the model are pre- 
sented. Furthermore, in order to complete this analysis, some energetics consideration are 
carried out in section 4. The last section summarizes our results and in it we make some 
suggestions for further research on the subject. 



II. DESCRIPTION OF THE MODEL 

To study the DNA torsional dynamics in the spirit of the model proposed by Yakushe- 
vich [^|6|,[7|,|8|,|9| Jl^ , pT]jr2] , p!3|] , we consider the DNA molecule as two long elastic and weakly 



interacting rods that are bound around each other, producing a double helix. Each rod 
represents one of the polynucleotide chains of the DNA molecule. For simplicity, we neglect 
the helical structure, and instead of the double helix we consider only parallel rods each 
having the form of the straight line |]5|J^J^,3,P, P^ipT| , P^iP^l • Assuming that the twist and the 
bend motions are independent and only taking into account the first of the above motions, 
we have the mentioned torsional model. 

One of the limitations of the model considered above is that in it the effects of inhomo- 
geneities are neglected. The changes in the strength of hydrogen bonds for different pairs 
of nitrogen bases, i.e., the inhomogeneous character of the DNA molecule, can be taken into 
account incorporating to the equations proposed by Yakushevich a new term that repre- 
sents the presence of such inhomogeneities. Assuming a local character for these bonds, 
distributed along the double helix in sites denoted by z n ,n = 1,2,.. the DNA torsional 
dynamics is now described by the following equations: 

h-^~ K ^~^k + k—[(2R 2 + l R) Bin(^) - R 2 sin^x + 2 )] 

+MJ2$(z ~ z n )^[{2R 2 + l R) Bin(^) - R 2 sin(0 1 + 2 )] = (la) 

h-^ - K 2 a 2 -^ + k—[(2R 2 + loR) sin(0 2 ) - R 2 anfa + cj> 2 )) 

+M S(z - z n ) ^[{2R 2 + l R) sin(0 2 ) - R 2 sin(0 1 + 2 )] = (lb) 



2 



where Al/l = 1 — Iq[(2R + Iq — -Rcos0i — _Rcos0 2 ) 2 + (-Rsin^i — -Rsin0 2 ) 2 ]~2 and 8(x) is 
the Dirac's Delta function of argument x. 

In these equations, <f>i(i = 1,2) is the rotational angle, I{ is the moment of inertia, Ki is the 
rigidity of the longitudinal springs of the i — th chain, k is the rigidity of the transversal 
springs connecting both chains(weak bonds), M is the rigidity of the springs modeling the 
local inhomogeneities that connect the chains in the positions labeled by z n , R the radius of 
the chains, Iq the maximum separation between chains and a is the characteristic length of 
the base pair in the double helix(see Fig.l). As in reference, we can simplify these equations 
by assuming Z = 0, which leads to: 

- K,a 2 ^ + LR^sin^) - sin(& + 2 )] 
+MR 2 8{z - [2 sin(0x) - sin(0! + 2 )] = (2a) 



d 2 



d t 2 d z 2 

+MR 2 K z ~ z n ) [2 sin(0 2 



+ kR 2 [2 sin(0 2 ) - sin(0i + <f> 2 )] 
sin(0! + 2 )] = 



(2b) 



III. ANALYSIS OF THE MODEL 

Now we proceed to analyze the Eqs.(^) in a general way. As it is of common usage. 
Let us to introduce the traveling wave variable £ = z — vt, where v is a constant. Then we 
obtain the following system: 

rj 2 / 

~ W * ajr + ^ 2 [2sin(^) - sin^ + cj> 2 )\ 

+MR 2 £ S(£ - £„) [2 sin(^) - sin^ + 2 )] = (3) 

n 

with i — 1, 2 

where Wi = Kid 2 — IiV 2 and £ n = z n — vt. In the last system, we have considered for 
simplicity the presence of one inhomogeneity located at the point z . 
The system Eqs.(3) is equivalent to the following dynamical system: 

$ (4a) 



|| = ^(kR 2 + MR 2 5(i - eo))[2sin(0 1 ) - sin^ + 2 )] (4b) 
* (4c) 



^ = ±-(kR 2 + MR 2 5(i - £o))[2sin(0 2 ) - sin^ + cj> 2 )\ (4d) 
For points £ ^ £o ; t ne Eqs.@ reduces to the system studied in reference fll7 |. 
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d(f) 2 



1 
1 



fci? 2 [2sin(0i) -sin(0i + 2 )] 



kR 2 [2 sin(0 2 ) -8111(0! + 2 )] 



The last equations corresponds to the potential function: 

V(4>i, 2 ) = kR 2 [2(cos(j)i + cos0 2 ) - cos(0i + 2 )], 



(5a) 
(5b) 
(5c) 
(5d) 

(6) 



which has a maxima at the point <p\ = 2nir, 2 = 2mir, with n = 0, ±1,±2, ... and 
m = 0, ±1, ±2, .... the points 0i = (2n — 1) 7r, 2 = (2m — 1) n, correspond to local minima, 
while the points <p\ = (2n — 1) 7r, 2 = 2m7r and 0i = 2ri7r,0 2 = (2m— 1) 7r are saddle 
points. The trajectories joining these critical points correspond to kink solitons |L7 |. For 
£ ~ £o the system has the same soliton solution, now perturbed by the presence of the local 
inhomogeneities at this point, as we will see later. 

Let us analyze the situation <pi = — 2 . In references [|17] , [T8| it has been shown that for this 
case there exists exact soliton-like solutions. If the distance between the inhomogeneity and 
the center of mass of the soliton is much greater than the soliton radius, the interaction 
between these two entities is very weak, and our model reproduces the solution of the 
mentioned references. The interaction becomes appreciable at a certain critical distance, 
that we will estimate later below. 

Let us consider symmetrical chains, i.e.;Ji = J 2 = / and K\ = K 2 = K and only one 
inhomogeneity. Then we can write Eqs.(@)as: 



1 d\ 



cf 2 



1 f^[2sin( 



d t 2 d z 2 A 2 

Q 



sin(0i + 4> 2 )] 



+ ^(z - zq) [2 sin^x) - sin^x + 2 )] = 



(7a) 



1 <9 2 2 
c 2 d t 2 



d 2 



jjr +^[2sin(0 2 )-sin(0 1 + 2 )] 

Q 



+ ^5(z-zo)[2sin(0 2 ) -sin(0 1 + 2 )] = (7b) 

where A 2 = Ka 2 /kR 2 , Q = M/k and c 2 = Ka 2 /I. 
To investigate the stability of non-perturbed soliton solution that in our case is [0,|T3], 



4>i(z) = — 4>2(z) = 4>k{.z) = 4arctan 



cxp 



Z Zr, 



2r 



(8) 



by 



where z m is the point that localizes the soliton radius r(half width of the soliton) defined 



A/2V2 



(9) 
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We consider small-amplitude excitations around <f>i(z), i.e., we put: 

(f>i(z, t) = <f>i(z) + fi(z) exp(Xit), i = l,2 



(10) 



where fi(z) are small. Linearizing in a standard way |19| Eqs.(0) around (pi(z,t) we get 
the following spectral problems for the exited modes: 



■0 + ^/i + ^[2oos(^)/ 1 -/ 1 -/ a ] 



+ ^6 (z - z ) [2 cos (<f> K ) h-h- A] = 



d 2 h , A 



dz 2 



+ ^/2 + ^[2cos(0 i ,)/ 2 -/ 1 -/ 2 ] 

+|^-* ) [2 cos (0k) / 2 - A -M 



(11a) 



(lib) 



We can write the above system as follows: 
d 2 F± 4 F + 



cfe 2 



N±F±+ 



A 2 cosh 2 (V2 (* - z m ) /A 
2 



2g 

A 2 



cosh 2 (\/2(z - z m ) / A 



5 (z - z ) = 



(12) 



where _F + and F_ are symmetric and antisymmetric solutions, i.e.,-F± = (fi ± A)/2, 
N + ,N_ corresponds to F + , F_ and: 



^2 ' 



A 2 2 



A 2 



(13) 



We can see that x = z o — z m represents the distance between the inhomogeneity and 
the soliton's center of mass. The soliton is unstable for positive eigenvalues. For the spec- 
tral problem Eqs. (|i~2f) , (|l3|)there exist two positive eigenvalues corresponding to the discrete 
spectrumfsee Appendix]. The solution corresponding to the first positive eigenvalue: 



A, 



1 2 



\ A 2 A 4 [cosh 2 ( v / 2x/A) 



(14) 



is 



exp 



A- 



Q 

A 2 " 



cosh 2 ( v / 2x/A) 



h = A = = 



exp 



A- 



Q 

A 2 " 



cosh(\/2(2; - z m )) 



cosh 2 (V2x/A) 



{z - zo) 



(z - Zq) 



cosh(y / 2(2; — z m )) 



Z > Z 



Z < Z Q 



(15) 
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where A is determined from the normalization condition. 
For the other positive eigenvalue: 



A 2 



cosh 2 (v^x/A 



- 1 



(16) 



we have for j\ = —f'2 = F-(z) the same expressions Eq. flT5|) . As we can see from 
the equations above if we put in these Q = and z m = 0, the results coincide with those 
reported in reference |T7J for these solitons. From the condition of instability (A > 0), the 
equation Eq. (|A2|) and the condition Eq. (|A6|) we find that for only distances smaller than 
a critical distance the inhomogeneity can affect the structural stability of the open state, 
i.e., that for distances(between the propagating soliton and the inhomogeneity) greater than 
such a critical distance, the interaction only influences the dynamics of its motion. This 
fact can be explained because for greater distances the inhomogeneity has not influence on 
the discrete states of the eigenvalue problem for the stability analysis, namely its solution 
reduces to the translational mode which is stable. This change in the stability can only to 
be due to a releasing of energy, i.e, a rupture of the stronger bond. 

The critical distance at which a rupture of the stronger bond occur, can be obtained from 
the Eq. (|A9|) . This can written as a function of soliton's radius as: 



X 



d, 



crit 



1.606r 



(17) 



The existence of a critical distance starting from which the bond is broken giving place 
to new implications no cleared up to now. For that reason we make in the following a simple 
energetic analysis. 



IV. SOME ENERGETICS CONSIDERATIONS 

In this section we briefly discuss some aspects concerning to the energetic contributions 
that appear in the presence of the inhomogeneities and it is a proposed field in order to 



model the motion of solitons as a point-particles as is usually done [f2(],[n],^,[24]] . The total 
energy for this mechanical model without to take into account the contribution due to 
inhomogeneities is given by the following expression: 

Eq = Ex + Eg + Et (18) 

where E g , Et are the energies due to elastic force between the cylinders of different 
chains and the torsional force between the cylinders of the same chain respectively. Ek is 
the corresponding kinetic energy. 
In the presence of inhomogeneities E has the form 

E = E + V (19) 

where the term 

V = MR2A4 > 2 (20) 
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takes into account the difference between the strength of weak and strong bonds. 
From the energetic point of view the traveling solitary wave must give an additional energy 
for breaking such a bond. Neglecting the dissipative effects and making the assumption that 
the rupture is produced when the solitary wave center is separated from the inhomogeneity 
by the critical distance then the rupture is possible if the effective kinetic energy has a value 
greater than the potential energy at the critical distance. 

If this hypothesis is true it is possible to model the motion of solitons as free particles in 
a potential field V with an effective kinetic energy E . This function V is a renormalized 
potential, which picking up the main effects due to the interaction with the inhomogeneity. 
In order to build the potential function we have taken into account that for distances lower 
than the critical distance between the center of soliton and the inhomogeneity the elastic 
energy becomes zero because of the rupture of the bond at this distance. In a general form 
the potential function is: 

Zq)) - - <fa.p{£ - ^o))] 2 i ~ Zq > d crit 

|£ — z o\ < d C rit 
Zo)) - (01u(-Oo) - 0l p (£ - Z ))] 2 £ - Z < -dcrit 

(21) 



v(0 





-^[(0 2u (-oc)-0 2p (£ 



The functions ip with % = 1,2 are the unperturbed functions Eq.(|S|), plus certain induced 
perturbations by the inhomogeneity. A suitable solutions for (p ip are possible to find by 



using the perturbative approach f2"2"|| . However, due to the extenuant calculus, we propose 
with good approximation the perturbation given by a function like the shape mode, first, for 
the relevant contribution that it plays in the interaction dynamics with the inhomogeneity, 



outlined by some authors in another generic problems [23] and second, because the main 
effect produced by inhomogeneity in its relative motion towards on solitons in this model is 
a contraction of soliton radius which is well represented by the shape mode, for being the 
inhomogeneity a simple bond between the chains, on the contrary of others similar models 
]2|1I0]'@]- Also we consider only the situation of a slow motion for the solitons which does 
not take into account the deformation due to the velocity. For the case of antisymmetric 
functions we get the following potential function: 



v(0 



2MB 2 


2MB 2 



2n — 4 arctan 



exp 



^ sinh(/^(C-Zo)) 
T cosh 2 (/3^%-^o)) 



4 arctan 



exp 



+ 



T cosh 2 (/^%-^ )) 



£ — ZO > d C rit 
|£ — z o\ < dcrit 
£ — z < ~ d cr it 

(22) 



where (3 = is a parameter that we consider gives a measure of the deformation due to 
the strength of the inhomogeneity. 
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In order to clarify our results, we will proceed to compare Eq.(|2lD with a potential function 
calculated by standard methods pOj^Tl . 

In our model the stronger bonds distributed along the chain correspond to static inho- 
mogeneities. The static inhomogeneities have been considered as perturbations in many 
nonlinear models [l6",pU,2"I|. By using a perturbative procedure we can estimate the effect 
of inhomogeneities on solitons dynamics if the characteristic size of the impurities is either 
much larger or much smaller than the soliton size. For take into account only the zeroth- 
order approximation we suppose that the influence of the inhomogeneities on the soliton 
gives rise only to a change in their parameters, like the center position, velocity, ect. It is 
usually named adiabatic approach 



20 



For the Eqs(0), taking the antisymmetric solutions, the coupled equations become: 



1 d 2 <p d 2 <p 2 

^dl 2 ~dl? + a 2 " 8111 ' 



2 ®U \ ■ 
~ A 2 " _ z o) sm ' 



(23) 



This is the perturbed Sine-Gordon equation which for Q = has a kink-like solution 
given by: 

Among the quantities that are conserved during the evolution there is the momentum 



00 d(j) dcj) 
-co d t d z 



(24) 



In the presence of inhomogeneities the momentum is no longer conserved by using 
Eq(4.6) it is possible to show that it varies according to |20| , 21| , |23 , |24 , 25] : 



dP 
~dt 



dz—6 (z 



d (cos < 
d z 



(25) 



The adiabatic theory is defined by the assumption that for small enough values Q, 
the kink shape will not be affected and its coordinate xi and its velocity v will become 
slowly changing functions of time. Within this hypothesis it can be shown that, in the 
non-relativistic limit v 2 c 2 , the kink center obeys the following evolution law. 



dt 2 ~ ' 

where the potential function is 

u(0 = 



dz—6(z 



d (cos < 
d z 



4Q 



A 2 cosh 2 (V2(£-Zo)/A) 



(26) 



(27) 



and £, zq are the positions of the mass center and the inhomogeneity respectively. 
The Eq. (P?D is a dimensionless function, multiplying it by the term Ka 2 this becomes 



u(0 



AMR 2 



COsh 2 (y/2 - Zq) /A 



(28) 



By comparing both Eq.(p2|) and Eq . (|28|) it is easy to(see Fig.2)that there exist a great 
differences between the profiles of the two curves. 
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Such differences are explained, on the one hand, because the calculation of the standard po- 
tential was done with the zeroth-order approach, which only produces changes in parameters 
like center position, velocity, without altering the shape. On the other hand, our proposal is 
a renormalized function which takes with a good approximation the deformation due to the 
interaction with the inhomogeneity. However the great difference owes its cause to the fact 
that the inhomogeneities on the contrary of other models, are localized in the bonds between 
the chains, therefore when the solitons (rupture of bonds in a localized region along the cou- 
pled chains) are situated from the inhomogeneity below the critical distance it is removed 
because the bond is broken. From the energetic analysis and taking into a account this last 
issue we derive two implications, first, the maxima elastic potential energy is reached for this 
critical distance and second, a potential well is formed whose width is given by twice the 
critical distance, such as is depicted in the (Fig. 2). The existence of this potential well could 
given place to interesting dynamics for example localized traps for the solitons distributed 
along the couple chains. Beside DNA being a low dimensional system these results could 
give many possibilities in its applications for the develop of new nanoelectronic devices as 
reported for such systems p6| , p7| . 



V. CONCLUDING REMARKS 



We have extended the model of Yakushevich [|13j for the DNA torsional dynamics, to 
take into account the presence of inhomogeneities in the molecular chains. The stability of 
certain kind of solitons that can propagate through the molecule, simulating the dynamics 
of open states, are studied. The obtained results are more general than those reported in 



the literature [|T7],|T3|] and include them as particular cases. We have found that for distances 
larger than 1.606 times the soliton's radius, the effect of the inhomogeneities in the structural 
stability of the open state can be neglected. Furthermore we analyzed the motion of DNA 
torsional solitons as the motion of point-particles through a renormalized potential function 
expressed in terms of the critical distance, which could give interesting dynamics. 
Our study has been done considering the presence of only one inhomogeneity and for the 
antisymmetric solitons reported in [HJ as exact analytical solutions for the unperturbed 
model. Subsequent studies can be directed to analyze the behavior of solitons between two 
inhomogeneities in terms of soliton's radius and the distance between inhomogeneities ratio. 
The more general case including a well-determined distribution of those bonds along the 
double helix will be considered in a future work. 
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APPENDIX 

The Eq.(12) is a eigenvalue problem which have a quantum mechanical analogue: 
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being: 



d 2 (f) 2m\ E\ 

dz 2 w~ 



fix) 



+ 



2b 2 



cosh 2 (y/2(z - z m ) /A 



+ 2Df( X )5(z-z ) = 



(Al) 



x = z - z n 



cosh 2 (V2x/A) 

where the potential function of the Schrodinger problem is: 



(A2) 



2b 2 



2m \ C osh 2 (y/2(z - z m )/A 



+ 2Df( X )S(z-z ) 



(A3) 



We propose the following solution: 



(*) = < 



' A exp\-Df(x)(z-z )] z>z 
cosh(6(z - z m )) 



A exp[Df(x)(z - zp)] 
cosh(6(z — z m )) 



(A4) 



z < z Q 



with the normalization condition 

exp[2Df(x)(z-z ) 



A' 



-oo cosh 2 (6(2; — z m )) 
whose eigenvalue is 



. dz+ , = 1 



cosh (b(z — z m )) 



(A5) 



This solution Eq.(Al) in the intervals (— oo, z — cr},[z Q + a, +oo) where a can take a very 
small value. If we put the proposed solution in Eq. flAlp , integrating in the interval and then 
tending to the limit a — > we can demonstrate the uniqueness of the solution with the 
constraint: 

fix) > (A6) 
For degenerated cases, i.e., D = 0, b ^ and D ^ 0, b = 0, the Eq. (|Al|) reduces to: 



d 2 d) 2m\ E\ 



and 



dz 2 h 2 



d 2 (b 2m\ E\ 



+ 



2b 2 



cosh 2 (V2(z - z m )/A) 



with E < 



(A7) 



dz 2 



2Df(x)S{z~ z ) = withE<0,D>0 



(A8) 



respectively 

The above equations were studied in reference |28|.It is easy to show that our proposed 
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solution Eg. (|A4|) with eigenvalues Eg. (|A5|) reduces to those proposed for the corresponding 
particular cases, i.e. Eg.(|A7|) and Eg. (|A8|) . 

By comparison of both Eg. ([12]) and Eg. ( |A1| ) using the above results we obtain the formulae 
Eg.©- Eg.©. 

For the case f(x) < 0, we would obtain a function potential with a similar shape but 
now with a reverse Delta function away from the well formed by the hyperbolic function 
potential. In this case the stationary states formed in the well of the potential function are 
not affected and corresponds with the eigenfunction for D = Eg. (|A7| ). 
From the condition 

fix) = fidcrit) = (A9) 

derives a distance starting from which either the well potential introduced by Delta function 
can have influence on the discrete states formed in the well of the function potential or not (see 
Fig3).The Fig.(3) is the plot of the potential function(Eg.(A3)) with z m = 0, b 2 = 1/2 and 
the rest of parameters that not appear in the figure are set to one for simplicity. 
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FIGURES 

FIG. 1. Section of coupled chains which depicts two couple strands of DNA with a stronger 
bond joining them. 

FIG. 2. Field Potential created by the inhomogeneity in which the soliton moves as 
a point-particle, a) Standard potential function (dotteddine), b)Renormalized potential func- 
tion(continuous-line) . 

FIG. 3. Potential function for a equation like-Schrodinger equation derived from the stability 
analysis. 
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